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Abstract 

In this article, we study canard solutions of the forced van der Pol equation in the relaxation limit 
for low-, intermediate-, and high-frequency periodic forcing. A central numerical observation made 
herein, which motivated our study, is that there are two branches of canards in parameter space which 
extend across all positive forcing frequencies. In the low-frequency forcing regime, we demonstrate 
the existence of the primary maximal canards induced by folded saddle-nodes of type I, and establish 
explicit formulas for the parameter values at which the primary maximal canards and their fold 
curves exist. Then, we turn to the intermediate- and high-frequency forcing regimes, and show that 
the forced van der Pol possesses torus canards instead. These torus canards consist of long segments 
near families of attracting and repelling limit cycles of the fast system, in alternation. We also 
derive explicit formulas for the parameter values at which the maximal torus canards and their fold 
curves exist. Primary maximal canards and maximal torus canards correspond geometrically to the 
situation in which the persistent manifolds near the family of attracting limit cycles coincide to all 
orders with the persistent manifolds that lie near the family of repelling limit cycles. The formulas 
derived for the folds of maximal canards in all three frequency regimes turn out to be representations 
of a single formula in the appropriate parameter regimes, and this unification confirms the central 
numerical observation that the folds of the maximal canards created in the low-frequency regime 
continue directly into the fold curves of the maximal torus canards that exist in the intermediate- 
and high-frequency regimes. In addition, we study the secondary canards induced by the folded 
singularities in the low-frequency regime and hnd that the fold curves of the secondary canards turn 
around in the intermediate-frequency regime, instead of continuing into the high-frequency regime 
as the primary maximal canards do. Also, we identify the mechanism responsible for this turning. 
Finally, we show that the forced van der Pol equation is a normal form type equation for a class of 
single-frequency periodically-driven slow/fast systems with two fast variables and one slow variable 
which possess a nondegenerate fold of limit cycles. The analytic techniques used herein rely on 
geometric desingularization, invariant manifold theory, Melnikov theory, and normal form methods. 

The numerical methods used herein were developed in lfT3lfT4l . 
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1 Introduction 

The forced van der Pol equation is a fundamental model for oscillatory processes in physics, electronics, 
biology, neurology, sociology, and economics. Possessing strong nonlinear damping effects, it is the 
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prototype of a forced relaxation oscillator, exhibiting slow and fast time scales, see ll5l 1^ [TOl 1^ 1^ [Ml 
[39l[40l[50l|43]l45l. The equations may be formulated as 

x' = y- fix), 

y'= ei—x + a + hcos9), (1) 

o' = oj, 

where the prime denotes the derivative with respect to the fast time variable r, fix) = ^x^ — x, and 
0 < e <C 1. The external signal, a + bcosO, models a time-periodic driving force with drive frequency 
w > 0. We refer to equation ([T]l as the fvdP equation. Throughout this article, we will work with the 
form of the system given by ([T]), as it allows us to explore the full range of forcing frequencies w > 0. 

A number of detailed studies of forced van der Pol equations have been carried out in the low- 
frequency forcing regime, u = 0(e), see |l5l|27l|30l|52l. We cite in particular the study [5] of a forced 
van der Pol system with low-frequency forcing, which presents a detailed analysis of the folded saddle 
singularities and their attendant canards. In the context of excitable systems (in particular, neuronal 
models), the folded saddle maximal canard plays the role of an excitability threshold manifold, locally 
dividing trajectories between those that jump at the fold to a different attracting manifold and those 
that do not ll4^ l57ll . This is also true in planar neuronal systems where solutions containing maximal 
canard segments correspond to excitability thresholds both in the case of type I neurons (integrators) and 
type II neurons (resonators) lIT^ ISTll . More generally, the canards induced by folded singularities, of 
folded node, folded saddle, and folded saddle-node types, have also been studied in models of neuronal 
dynamics ||47]|48l|53l and in many other systems, see for example lfflr7l [T3l[T5ll37ll54ll55ll5^l57l . 

In this article, we examine the fvdP equation ([T]l in three different regimes of forcing frequencies: 
low-frequency (w = 0(e)), intermediate-frequency (tu = Oi^/e)), and high-frequency (io = 0(1)). In 
each regime, we study the canard solutions that the fvdP equation ([T]l exhibits. 

We begin in the low-frequency regime. First, we briefly apply the theory of folded singularities to 
([T]), to identify the different types of folded singularities that it exhibits in this regime. We place special 
emphasis on folded saddle-node singularities of type I (FSN I), which are known to generate a number 
of different types of canard solutions. 

The graph of the fast null-cline, y = fix), of system Q with e = 0 plays a central role in under¬ 
standing the system dynamics. We are especially interested in the repelling branch in the middle and the 
attracting branch on the right, respectively, of the graph. Let Sr denote the two-dimensional manifold 
formed by rotating the (middle) repelling branch through one complete revolution in the angle 9 : [0, 27r), 
and similarly let Sa be the the two-dimensional manifold formed by rotating the (right) attracting branch 
through one complete revolution in 9. In the low-frequency regime of ([T]l, Fenichel theory ll2^l^ guar¬ 
antees that, when e is sufficiently small, there exist two-dimensional, locally-invariant manifolds and 
near Sr and Sa, respectively, away from the fold regions. In the low-frequency forcing regime of ([T]), 
these persistent manifolds are referred to as slow manifolds, since the dynamics on them is slow in y and 
9. 

The primary canards of folded singularities are orbits that have a long segment close to S%, pass 
through a neighborhood of the folded singularity, and then have a long segment near S^. The lengths 
of these segments depend on the parameter values; and, there are curves of parameter values along 
which the segment near Sr has maximal length, going all the way up to the other fold curve. These 
primary canards are referred to as primary maximal canards. See Figure[TJa) for a representative primary 
maximal canard. They are determined geometrically by the parameter values for which Sr and agree 
to all orders in e, in a manner that is analogous to the maximal limit cycle canards (also known as the 
maximal headless ducks) in the classical, planar van der Pol equation, recall li^ fTTlI^ . 

The following is the first main result of this article: 

Theorem 1.1. low-frequency forcing. Let oo = eoJ, where a; = 0(1), and let b = O(y^). Then, there 
exists an eo > 0 such that for all 0 < e < eo, there are two curves in the (a, w) parameter plane given by 

a = 1 - I ± 6exp , (2) 
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emanating from the points (w, a) = (0,1 — | ± b), along which the system ([T]l has folds of primary 
maximal canards. Moreover, for each 0(1) value of uJ, the system o has two primary maximal canards 
for every value of a in the interval between the points on these fold curves. Finally, there are no primary 
maximal canards for values of a outside the closures of these intervals. 

This first theorem is established by using the geometric desingularization method, also known as the 
blow-up method |[T^ l20l 1^ . to inflate the folded saddle node points of type I into hyperspheres, and 
then by employing invariant manifold theory and Melnikov theory in the appropriate coordinate charts. 

Next, we show that the fvdP system ([T]l has torus canards in both the regime of intermediate- 
frequency forcing, in which the system ([T]l has three time scales: x is a fast variable, 6 is an intermediate 
time variable, and y is a slow variable; and in the regime of high-frequency forcing in which Q is a two- 
fast (x, 6) and one-slow system (y). Torus canards are a relatively new type of canard solution discovered 
in a model of neuronal activity in Purkinje cells O^ . They consist of long segments near families of 
attracting and repelling limit cycles of the fast system, in alternation. Torus canards have recently been 
shown to exist in a broad array of models, including in three models of neuronal bursting, see [81: the 
Hindmarsh-Rose model, the Morris-Lecar-Terman system, and the Wilson-Cowan-Izhikevich equations; 
in a model of elliptic bursters, where the torus canards are rotated versions of limit cycle canards of a 
planar system, see ll^ ; in a rotated van der Pol-type model system, see ||2l ; as well as more recently 
in a model of respiratory rhythm generation in a pre-Botzinger complex, see ll46ll . The significance of 
torus canards is that they play a central role in the transition between periodic spiking and bursting, of 
different types, in these neuronal models. 

(a) (b) 




Figure 1: (a) Segment of a primary maximal canard solution of Q, and (b) segment of a maximal torus canard 
solution of Q. Both have long segments near the family of attracting limit cycles (outer portion of the green 
surface) and near the family of repelling limit cycles (inner portion of the green surface). Here, a = 0.997, b = 
0.994, e = 0.02, and (a) uj = 0.001 and (b) w = 0.3. 

In the intermediate- and high-frequency regimes of ([T]l, Fenichel theory ll2^l^ also guarantees that, 
when e is sufficiently small, there exist two-dimensional, locally-invariant manifolds near Sr and Sa, 
away from the fold regions. We again denote these by and and label them as persistent manifolds. 
However, it is crucial to observe that these persistent manifolds are no longer slow manifolds in these 
regions. Instead, the orbits of ([T]l on these persistent manifolds exhibit two time scales, with fast rotation 
due to the limit cycles, as well as slow drift in the vertical direction, down and up along S^. 

Torus canards are orbits of ([T]l in the intermediate- and high-frequency regimes that have long seg¬ 
ments near S^, spiral through a neighborhood of the fold curve of limit cycles, and then have a long 
segment near S^- The lengths of time that torus canards spiral around near and depend on the 
system parameters, and for system ([T]l there are curves of parameter values along which the time spent 
near is maximal, with the orbits spiralling all the way up S^- These are defined to be maximal torus 
canards, in analogy with the maximal limit cycle canards of the unforced van der Pol oscillator. A 
representative maximal torus canard is shown in Figure [T](b). 

For system ([T]) in the intermediate-frequency regime, we prove the following theorem: 
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Theorem 1.2. intermediate-frequency forcing. Let w = ^/e^, where ri = 0(1), and let b = 0{e). 
Then, there exists an eo > 0 such that for all 0 < e < Eq, there are two curves in the (a, fl) parameter 
plane given by 


a = 1 — - ± 6 exp 
8 



(3) 


along which the system ([T]l has folds of maximal canards. Moreover, for each fixed 0(1) value of Q, 
the system ([T]l has two canards for every value of a in the interval between these fold curves, and none 
outside the closure of these intervals. 


This theorem is also established by using the geometric desingularization method; however, in this 
regime, we inflate the circular fold curve along which the attracting and repelling limit cycles meet into a 
two-torus, rather than the FSNI points. Also, the scalings are different, as is the analysis in the coordinate 
charts near the torus. 

Then, for the high-frequency regime, we establish: 

Corollary 1.3. high-frequency forcing. Let lo = 0(1), and let b = 0(e). Then, there exists an sq > 0 
such that for all 0 < e < eo, there are two curves in the (a, uj) parameter plane given by 


a = 1 — - ±b exp 
8 



(4) 


along which the system ([T]l has folds of torus canards. Moreover, there exists a pair of torus canards for 
each parameter value in the interval between these fold curves. 


We note that the presence of the torus canards in this type of fast-slow system is signalled by the 
existence of a fold of limit cycles of the fast system, here at (x, y) = (1, —|), together with a nearby 
torus bifurcation in the full system, here at 


1 — a 


2 


1 6 ^ 

2 (a^ — l)2ti;2 + (e —uj‘^)‘^ 


(5) 


See Appendix A. These two triggering mechanisms arise ubiquitously in fast-slow systems with two or 
more fast variables. 

Having established these theorems for the existence of the primary maximal canards and the torus 
canards, as well as their folds, we now analyze the relationship between these results. Plainly, the 
formulas for the curves of folds Q, Q, and Q in the three different regimes are all representations of 
the same formula. 


a = 1 — - ±b exp 



( 6 ) 


in the respective frequency regimes. The exponential term has magnitude b (which is 0(^/s)) and is 
slowly varying with u] in the low-frequency regime (Theorem [nj, small-amplitude (b = 0(e)) and 
varying with 0(1) frequency in the intermediate-frequency regime (Theorem |L2[ ), and exponentially 
small in e in the high-frequency regime (Corollary |L3| ). 

The analysis in all three regions shows that the values of the parameter a for which the canards exist 
in between the fold curves may similarly be summarized succinctly in one formula: 



6 cos( 00 )exp 



(V) 


Here, 6q is an arbitrary phase, and the magnitude and dynamics of the exponential term are also as 
discussed above. 

It is also of interest to observe that, in the limit of a; —)■ cx), formulas Q and (|7]l become a —)• Oc := 
1 — |, which corresponds exactly to the leading order locations of the maximal limit cycle canards in 
the planar vdP equation, see for example 1361. Therefore, as expected, for sufficiently high-frequency 
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Figure 2: Curves of folds of maximal canards in the (w, a) plane as obtained from Q (red curves) and 
2-parameter numerical continuation (blue curves) for e = 0.01 and (a) b = 0.01; (b) b = 0.02; (c) 
b = 0.0351 and, (d) b = 0.1. For b = 0{e) (panels (a)-(c)), there is good agreement between theoretical 
(red) and numerical (blue) results over the entire range of forcing frequencies, including for both the 
primary maximal canards which exist for io = 0(e) and the maximal torus canards which exist for 
u! = 0(1). Note that the scales in panels (a)-(c) are the same. For b = 0{^/e) (panel (d), in which the 
vertical scale is different), we find that the numerical continuation terminates when ui is no longer 0(e). 


forcing, the effect of the forcing averages out to this order, and ([T]) behaves like the classical planar vdP 
equation. In this limiting regime, the torus canards of ([T]l appear to be rotated copies of the limit cycle 
canards of the planar vdP. 

Then, with the above analytical results in hand, we turn next to the results of numerical continua¬ 
tions which confirm that the curves of the folds of primary strong canards observed in the low-frequency 
regime continue directly to the curves of the folds of torus canards discovered in the intermediate- and 
high-frequency regimes. This is illustrated in Figure Moreover, as is also shown in this figure, the 
agreement between the formulas and the numerical continuation results are excellent within the parame¬ 
ter regions stated in the theorems. We also note that the theory does not appear to extend outside of these 
regions, and also preliminary results of numerical continuation reveal different dynamics there. Overall, 
then, the formulas Q, and Q together with the numerical continuations, will directly imply that the pri¬ 
mary strong canards, which exist in the low-frequency forcing region, continue naturally to the branches 
of torus canards, which exist in the high-frequency regime, where the folded singularities cease to exist, 
with the transition happening in the intermediate-frequency regime. Understanding the continuation dy¬ 
namics of these curves is one of the main results of this article. Moreover, the results here will also help 
shed light on other models with torus canards. In particular, we observe that numerical simulations of a 
rotated van der Pol-type model exhibit the same continuation of the maximal canards across the entire 
range of forcing frequencies; see Figure 5 in [2]. Numerical continuations in other neuronal (or neural) 
models lUl show similar phenomena. 

In this article, we also study the secondary canards of ([T]l. Secondary canards, lie near the primary 
canards for most of their lengths, and make a number of small loops around an axis of rotation usually 
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Figure 3: Curves of folds along the computed branches of primary and secondary canards. Here, e = 
0.05 and b = 0.01. 


referred to as the weak canard. Secondary canards are indexed by the number of loops they make around 
the weak canard and by the value of the y-intercepts of the solutions during the nearly-horizontal jumps 
that occur from a neighborhood of the family of repelling limit cycles back to a neighborhood of the 
family of attracting limit cycles. In particular, we carry out numerical continuation of the branches of 
maximal secondary canards of ([T]) that are created by the folded singularities in the low-frequency regime. 
In contrast to what we find for the primary canards, the branches of the secondary canards turn around 
well before they reach the high-frequency regime. See Figurej^ Also, we identify the mechanisms which 
cause the branches to turn. 

To conclude this article, we demonstrate that ([T]l serves as a local normal form for slow/fast systems 
with one slow variable and two fast variables in which the fast subsystem possesses a non-degenerate 
fold of limit cycles, and in which the slow system is subject to time-periodic forcing. These fast-slow 
systems exhibit torus canard explosions, just as shown here for the fvdP equation ([T]l, and one may 
therefore directly conclude, by applying the same techniques used herein, that the folds of their canards 
behave in a similar fashion. 

The numerical method developed in ifT^ [T4ll . and also employed in ifTSl . is the main numerical 
method used throughout this article to find the persistent invariant manifolds and the curves of maximal 
canards that lie in the intersections of these manifolds. This method, which uses the AUTO continuation 
software m, turns the problem of finding the invariant manifolds of fast-slow systems into a boundary 
value problem for system ([T]l with the time of integration included as a parameter. Then, the parametrized 
families of solutions of the two-point boundary value problems are continued. This allows to integrate 
in positive and negative time using pseudo-arclength continuation, approximating the orbit segments of 
solutions of system ([T]) subject to particular boundary conditions by orthogonal collocation, which is 
very well suited to multiple time scale vector fields (see MM)- 

The outline of the article is as follows. In Section we consider system ([T]l with low-frequency 
forcing, cc = 0{e), and apply the theory of folded singularities in a straightforward manner to find 
fhe associated primary and secondary canards of folded singularities. Then, in Section we prove 
Theorem o establishing the existence of the primary maximal canards induced by the FSN I points, 
and the associated fold curves, including the derivations of formulas Q and ([7]) for system Q in the low- 
frequency 0{e) region. We then turn to the cases of intermediate- and high-frequency forcing in Section 
where we study the torus canards of Q. We prove Theorem 1.2 and Corollary |1.3[ establishing the 
existence of the maximal torus canards and their fold curves in these regimes, as given by the formulas Q 
and Q. This shows analytically that the curves of the fold curves of the primary maximal canards, which 
are bom in the low-frequency regime, continue for all a; > 0 into the fold curves of the maximal torus 
canards. We also observe that the analytically-derived formulas and the curves obtained in the numerical 
continuations agree over the entire range of forcing frequencies. Then, in Section we examine the 
numerical continuations of the folds of secondary canards, and we identify the mechanism by which 
they turn around well before they reach the high-frequency regime. Also, we analyse how the curves 
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Figure 4: 2-parameter continuation of folds of primary maximal canards of a folded saddle-node (type 
I) for e = 0.01 and b = 0.01. Orbit segments are plotted in ‘Cartesian’ coordinates {u, v, y) = 
{xcos9,xsm6,y). 


of the folds of secondary canards induced by folded nodes change as the parameter h is varied, up to 
and including b = 0(1). and hence as the distance between the folded node and the folded saddle is 
varied. The final main result of this article is presented in Section We demonstrate that o may be 
considered as a local normal form for some generic fast-slow systems that have a fold of limit cycles 
and that undergo a torus canard explosion. In Appendix A, we prove, using second-order averaging, the 
existence of a torus bifurcation in Q and calculate the locus Q for this bifurcation in parameter space. 

2 Low-Frequency Forcing: Canards of Folded Singularities, Especially 
of FSNI Points 

In this section, we present a brief review and analysis of the folded singularities that system ([T]l possesses 
in the regime of low drive frequency, i.e. uj = slJ, where uJ = 0(1). Readers familiar with the theory of 
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folded singularities and their canards may proceed to Sectionj^ In this regime, ([T]) is 

X = y- f{x), 

y = e {—X + 0 + 6 cos 9) , (8) 

9' = eLJ. 

It is a l-fast/2-slow problem with fast variable x and slow variables {y,9). We analyse the reduced 
dynamics associated to @ and derive the desingularised vector field on the critical manifold. Then, we 
identify the canards of the folded singularities. 

2.1 The Layer and Reduced Systems 

Taking the singular limit e —)■ 0 in ([^, one finds fhe ID layer problem 

x' = y-f{x), (9) 

where y and 9 are parameters. Alfemafively, fhe singular limif e —)■ 0 in Q gives fhe 2D reduced system 

0 = y- f{x), 

y = —X + a + bcos9, (10) 

9 = oJ, 

where the overdot denotes the derivative with respect to the slow time t = er. The manifolds 5^ and 
are non-unique. Hence, the canards that lie near the manifolds are also non-unique. However, for a fixed 
choice of invarianf manifolds, and fheir fransverse infersecfions correspond fo maximal canards. 

Sysfems Q and ( fTO] ) are fwo differenl approximations of fhe forced vdP equafion. The idea of 
Geomefric Singular Perfurbafion Theory (GSPT) ll^ is fo combine informalion from fhe ID layer and 
2D reduced problems in order fo undersfand fhe dynamics of fhe full 3D forced vdP equafion for 0 < 
e <C 1. 

We begin wifh an analysis of fhe ID layer problem Q, which is an approximafion of ([T]) wherein 
fhe slow processes are assumed fo move so slowly fhaf fhey are fixed. The critical manifold is fhe sef of 
equilibria of fhe layer problem ([^: 

S ■■= {{x,y,9) e X ; y = f{x)} . 

Linear sfabilify analysis of Q shows fhaf fhere are disjoinl curves, L, of fold poinfs given by 

L := {{x,y,9) e S : x = ±1} , 

which separafe fhe outer affracfing sheefs, Sa, of S from fhe middle repelling sheef, Sr, of S. Fenichel 
fheory ll2^ |33l guaranfees fhaf fhe normally hyperbolic segmenfs of S (i.e. fhe parts of Sa and Sr af 
0(1) disfances from fhe fold curve L) will persisf as invarianf slow manifolds, S^ and S^, of ([T]l for 

0 < e <C 1. 

The price we pay for fhe approximation Q is fhaf we have frivial dynamics on S. To obfain a non- 
frivial flow on fhe crifical manifold, we fum fo fhe reduced problem. The 2D reduced problem ( fTO] ) is an 
approximafion of ([T]l wherein fhe fasf motions are assumed fo be so rapid fhaf fhey immediately settle fo 
fheir sfeady sfafe, which is precisely fhe crifical manifold. In ofher words, fhe reduced problem prescribes 
a non-frivial flow along S. The price we pay for fhis particular approximation is fhaf fhe reduced flow is 
nol defined away from fhe crifical manifold. Note fhaf fhe resfricfion of fhe flow of ([T]) fo S^ is a small 
smoofh perfurbafion of fhe reduced flow on S. 

To analyse fhe flow on a manifold, we use fhe coordinafes {x,9). We differenfiafe fhe algebraic 
consfrainf y = f(x) wifh respecf fo t fo obfain fhe evolufion equafions in fhis coordinate chart, 

— l)i; = —X + 0 + 6 cos 9, 


9 = oj. 


( 11 ) 


The reduced flow (dT) is singular along the fold points L of Q. To remove the finite time hlow-up of 
solutions at the folds, we rescale time dt = {x^ — 1) ds to obtain 


X = —X + a + b cos 9, 

6 = co(x^ — 1 ), 


( 12 ) 


where the overdot now denotes derivatives with respect to s. System ( [T2l ) is equivalent to the reduced 
flow on the attracting sheets Sa, where the time rescaling dt = — l)ds preserves the orientation 

of trajectories. On the repelling sheet Sr, however, we have — 1 < 0, so that the time rescaling reverses 
the orientation of trajectories. Thus, to obtain the reduced flow ( [TT] ) on Sr from ( [T^ , we simply reverse 
the direction of trajectories of whenever we are on the repelling sheet of the critical manifold. 


2.2 Folded Singularities & Singular Canards 

The desingularised system ( [T2l ) possesses special equilibria called folded singularities, M, which are 
points along the fold curves where the right hand side of the x-equation vanishes. In system ( [T2l ), there 
are infinitely many pairs of such points (when 6 is considered in its lift to M): 

M := I (x, y, 0) G L : 0 = Ikir ± cos“^ ^^ , /c G Z 

where |1 — a| < b. Folded singularities are not true equilibria of the reduced flow ([TT]). Instead, they 
correspond to points of ( [TT] ) where there is potentially a cancellation of a simple zero in the x-equation 
and trajectories may pass through the fold (via the folded singularity) with finite speed. Such a trajec¬ 
tory of the reduced flow that passes through a folded singularity and crosses from the attracting (resp. 
repelling) sheet to the repelling (resp. attracting) sheet is called a singular canard (resp. singular faux 
canard) ll5Tll55ll5^ . 

Considered as equilibria of the desingularised system ( [T^ , folded singularities are classified accord¬ 
ing to their linearization. Folded nodes have real eigenvalues of the same sign. Folded saddles have real 
eigenvalues of opposite sign, whilst folded foci have complex eigenvalues. In the forced vdP equation 
Q, we find that for cu > 0 the folded singularities with 

6s{k) = 2k-K — cos“^ 

are folded saddles, whilst the folded singularities with 

9n{k) = 2kTT + cos“^ 

are folded nodes provided 

[l-a)^ <b^ + 

Folded nodes and folded saddles have been demonstrated to be the organising centres for complex phe¬ 
nomena. Folded nodes for instance have been identified as the cause of the small oscillations in MMOs 
in various neurophysiological problems |[22l . such as in a self-coupled FitzHugh-Nagumo model, in a 
Hodgkin-Huxley model |f48l and in a pituitary lactotroph cell model Il5^ . More recently, folded sad¬ 
dles have been identified as playing a significant role in distinguishing between transient spiking and 
quiescence in a model of propofol anaesthesia Il42]l and in non-autonomous excitability models ll57ll . 





2.3 Canards of Folded Saddle-Node Type I Points 

Folded nodes and folded saddles can be created through bifurcations in two distinct ways in ([T]): via a 
folded saddle-node (FSN) of type I |[54l[J7l or via an FSN of type II 071 . Both FSN’s correspond to a 
zero eigenvalue of the folded node (or folded saddle). The two FSN scenarios are distinguished by their 
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geometry. In the FSN I limit, the center manifold of the FSN I (in system ( [T^ ) is tangential to the fold 
curve. In the FSN II limit, the center manifold of the FSN II point is transverse to the fold curve. We 
focus here on FSN I and refer to the remark below for FSN II points. 

The FSN I is the codimension-1 bifurcation of the desingularised system ( [T^ in which a folded 
saddle and a folded node coalesce and annihilate each other in a saddle-node bifurcation of folded sin¬ 
gularities. For the forced vdP equation ([T]l, there are infinitely many such FSN I points: (x, y, 6) = 
(1, — 1,2A;7r), and they occur for a = 1 ± 5 and uj = ecJ. The FSN I at a = 1 — 6 has its center manifold 
on Sr so that the funnel region (enclosed by the strong canard of the folded node / folded saddle canard 
and the fold curve) vanishes in the FSN I limit. In this case, we expect generic solutions of 0 near this 
FSN I limit to either be relaxation oscillations or MMOs. The FSN I at a = 1 -|- 6 on the other hand 
has its center manifold on 5+ so that the funnel persists in the FSN I limit and most solutions of Q can 
tunnel through Sr and return to Sa- A representative example of the passage through a FSN I bifurcation 
at a = 1 -I- 6 is shown in Figure]^ 



Figure 5: Reduced flow of the fvdP equation a shown in a neighbourhood of the upper fold curve L (defined 
by X = 1, y = — |) for uj = l,b = 0.01, and (a) a = 1 -I- (b) a = 1 + 6 and (c) a = 1 -f In panel (a) where 

a < 1 -I- 6 , there is a folded node (0„) and a folded saddle (6s). The strong and weak eigendirections of the folded 
node are denoted by 7 ^ and 7 ^,, respectively. Similarly, the singular canard and faux canard of the folded saddle 
are labelled 7 s and 7/, respectively. Note that there is a heteroclinic connection CM from to 0s- In particular, 
7 u, is tangent to CM at 0 „, and 7/ is tangent to CM at 0 s- In panel (b) where a = 1 + b, the folded singularities 
merge to a FSN I (indicated by 0fsn)- In this case, the singular strong canard of the folded node merges with the 
maximal canard of the folded saddle. Meanwhile, the singular weak canard of the folded node merges with the 
faux canard of the folded saddle to become the center manifold of the FSN I. In panel (c) where a > 1 + b, 
the folded singularities (and associated singular canards) have been destroyed in the FSN I bifurcation, and there 
are no canard dynamics. 

In the FSN limit, the standard folded node/folded saddle theory requires modification. For the FSN 
I, the following results have recently been proved in IBTI . valid for 0 < e <C 1 and /r = 0(e") where 
a > 1/4: 

1. The singular strong canard of the folded node perturbs to the primary maximal strong canard. The 
singular canard of the folded saddle perturbs to a maximal canard. 

2. There exists a heteroclinic connection CM between the folded nodes and folded saddles of ( [T^ . 
This heteroclinic perturbs to a canard-faux canard solution CM^ that corresponds to both the 
primary weak canard of the folded node and the faux canard of the folded saddle (faux canards are 
the equivalent of singular faux canards for 0 < e <C 1). 

3. There exist canards and faux canards. 

Thus, canards and faux canards of the FSN I oscillate about an axis of rotation, which is approximately 
given by the heteroclinic CM that connects the folded node and folded saddle (see Figure |^a) for 
instance). For the fvdP, we find that CM := {(x, y, 0) G S' : x = a -|- 5 cos 9} . We study the associated 
maximal canards of 0 in Section]^ 
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Remark 1. For the forced vdP equation ([T]), FSN II points are codimension-2 bifurcation points of the 
desingularised flow (corresponding to tJ = 0), and constitute a special case of the FSN I. They can be 
analyzed using the approach presented in liTTl . 


3 Loci of the Maximal Canards for Low Forcing Frequencies 

In this section, we analyse system ([T]) with low-frequency forcing (w = suJ). We prove Theorem |1.1[ 
demonstrating that, for b = O^y/s), formula Q gives the branches of the folds of the primary maximal 
canards and that for each value of a in between the fold curves, there are two primary maximal canards 
of system o given in parameter space by ([T]). More precisely, we analyze the FSN I points to show that 
formula Q gives the locus of points at which the primary strong canard of a folded node point merges 
with the folded saddle maximal canard. We present the analysis for the FSN I that occurs for a = 1 — 6 
and note that the FSN I at a = 1 -|- 6 is treated similarly. 

For the analysis with low-frequency forcing, we first translate the FSN I at a = 1 — 6 to the origin 

2 

u = x-l, V = y +-, r] = a-l + b, 

O 

so that o is transformed to 

u' = v - -h , 

v'= £ (—u + rj + b(cos 0 — 1)), 

6' = £Uj. 


We then inflate the FSN I singularity to a hypersphere using the spherical blow-up transformation: 

u = r^M, V = f^v, 6 = f0, e = r^e. 


Moreover, we rescale the parameters b and 77 as 

b = y/£l3, 7? = £7, 

where /3 = 0(1) and 7 = 0(1). Also, we append the trivial equation e' = 0 to system ( [T^ and take 
77 > 0 sufficiently small. The (spherical) blow-up transformation is a map from B := x [—y, y] into 
M^. We examine the vector fields induced by this coordinate transformation in two useful coordinate 
charts: the entry-exit chart (or phase-directional chart) Ki : {v = 1} and the rescaling (or central) chart 
K 2 : {e = 1 }, beginning with Ki. 

In chart Ki, the blow-up coordinates are 

u = r\ui, V = rf, 0 = ri0i, £ = rf£i, (14) 


where the subscript corresponds to the chart number. The governing equations are 


2 1 2 s 1 

ill = 1 - Ui - 3^i^i ~ 2 ^^ 

ri = £i F, 

01 =ri£iuj -^£1 0iF, 

£1 = -£i F, 


(15) 


where F{ui, 0i,ri) = —ui + r‘l£i^ + f^^/Fi (cos(ri 0 i) — 1 ), and we have desingularised the vector 
field by a factor of rf and recycled the overdot to denote the derivative with respect to the new time. The 
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hyperplanes {ri = 0} and {ei = 0} are invariant. In the invariant subspace {ri = 0}, = 0 is an 

attracting fixed point. The line 


4 = {(ui,ri, 6 »i,ei) = (ui, 0 , 0 , 0 )} 

is invariant, and on it the system dynamics are governed by tti = 1 — Furthermore, on £u there are 
attracting and repelling fixed poinfs pa = ( 1 , 0 , 0 , 0 ) and pr = (— 1 , 0 , 0 , 0 ), which, respecfively, have 
cenfer manifolds A^a,i and Nr^i in fhe half space £1 > 0 . 

In order fo demonsfrafe fhe exisfence of fhe primary maximal canards, we will show fhaf fhere is a 
heferoclinic connecfion befween pa and pr in fhe hyperplane {ri = 0 } and fhaf Ibis heferoclinic orbif 
persisfs for sufficienfly small values of ri, using Melnikov fheory. The persisfenf connecfions correspond 
fo fhe primary maximal canards. We carry ouf fhe relevanf analysis in fhe rescaling charl K 2 , where fhe 
blow-up fransformafion is given by 


u = r2U2, V = r2V2, 0 = r202, e = ( 16 ) 

Nofe fhaf r 2 = so fhaf chart K 2 corresponds fo an e-dependenf rescaling of fhe forced vdP equafion. 
Also, fhe coordinafes in fhe fwo charfs are relafed via fhe following fransformafion: 

r 2 — nei , U2 — U1S1 , V2 — , ^2 — 016-^ , 


where ei > 0 . 

In charf K 2 , the blown-up system ( ffo] ) is 

U2=V2-ul- 

V2 =-U2 + rl'y + (3 {cos{r202) - 1), 

62 = r 2 w, 

where we have desingularised the vector field (i.e., rescaled by r^) and again recycled fhe overdof, fhis 
lime fo denole fhe derivafive wilh respecf fo fhe new lime t 2 - This syslem is singularly perlurbed wilh 
fasl variables {u 2 ,V 2 ) and slow variable 62 . Rewriting fhe blown-up syslem in non-aulonomous form, 
we have 


U2 = V2 -ul- ^r2ul, 

V2 = -U2 + r^-y -I- 13 [cos{r\u]t2) cos(r202,o) - l) - sin(r|wf2) sin(r26*2,o), 


(18) 


where 4,0 is an arbilrary phase. This is fhe syslem we will analyse in chart K 2 . For uj and t 2 of 0(1), 
we have 


cos{rlcjt2) cos(r24,o) - 1 = £>(^’ 2 ) ^s £2 —^ 0, 

sin(r 2 wf 2 ) sin(r24,o) = Oir^) as r 2 0. 


The unperlurbed problem corresponding fo ( fT 8 ) is oblained by selling r 2 = 0, 

u '2 = V 2 - ul, 

v '2 = -U2. 


( 19 ) 


This syslem is Hamillonian wilh Hamiltonian function 

H(u2, V2) = ^^2 - U2 - 0 , 
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Figure 6: Contours of the Hamiltonian function H. Periodic solutions of ( [T^ correspond to level sets 
with — H < 0. Unbounded solutions of ( [79] ) correspond to H > 0. The H = 0 contour, F, is the 
heteroclinic connecting the points pa and pr at infinity. 


and non-canonical formulation 


U 2 


V 2 


1 2 ..^ 

2 dv2' 


2 


dH 
du 2 ' 


The level curves of H are presented in Figure The contour F separates closed trajectories from 
unbounded orbits and has the explicit time parametrization 


(u2,r, ^^2,r) 



The separatrix F corresponds to the singular strong canard of the FSN I. In geometric terms, it is a 
heteroclinic orbit that lies on the upper-hemisphere and connects the fixed poinfs pa and pr, which bofh 
lie on fhe equafor of fhe blown-up sphere. 

We now use fhe Melnikov mefhod fo analyse fhe persistence of F under small-amplifude perfurba- 
fions. As applied fo ( fT^ , Melnikov fheory measures fhe splitting disfance D befween fhe curves of 
solutions of fhe perfurbed system fhaf are forward and backward asympfofic fo pr and pa, respecfively. 
We develop D in an asympfofic series in fhe small paramefer r 2 : 


D{r2) = dirl + d2rl H-, 


where fhe ferms in fhe Melnikov infegral are given by 


= 


V H\ 


_1,,3 

3 ^2 r 


^ 17 + ^ (cos(r|wf2) cos(r 26 » 2 ,o) “ l) j 


d 2 = 


fOO ( 0 \ 

y_oo ^ • I sin(r|a;f2) sin(r 26 ' 2 ,o) ) 


We note fhaf fhe infegrand in d 2 is an odd function of t 2 so fhe infegral evaluafes fo zero and fhe sine ferm 

has no confribufion fo fhe disfance measuremenf D. We also nofe fhaf fhe cos(r 2 e 2 ,o) i 

U 

di is 0(1) wifh respecf fo r 2 . The infegral di was evaluated by faking cos z = Re(e*^), completing fhe 
square in fhe exponential, and deforming fhe contour in fhe complex plane. The resulf is 


di 



P-rl 



(de 


2 2 ^ 


cos(r 26 ' 2 ,o) 
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Substituting this into the bifurcation equation D = 0, we have 


^2 


^+7 


+/3cos(r26'2,o) (^e 2 ’ 


- 1 = 0 . 


Thus, reverting to the parameters a, b, e, and to, we see that the primary maximal canards for the FSN 
I are given by 


a = 1 — - — bcos(0o) exp ( — 




( 20 ) 


which is Q. We remark that 0q is the arbitrary phase 02,0 in the original 0 coordinate (i.e., 0o = r 202 ,o)- 
This completes the demonstration that Q holds for low-frequency forcing, giving the locus of points 
at which the primary maximal canards exist. Moreover, one also sees that, for each uJ, the folds of the 
primary maximal canards at the endpoints of these parameter intervals are given by 


a = 1 — 


e 

8 


± 6 exp 



( 21 ) 


which is precisely 0 - The loci in the (w, o) plane of the folds of maximal canards mark the upper and 
lower boundaries of the regime in which the primary maximal canards exist. This completes the proof of 
Theorem ll.il 

Remark 2. For w = e tJ, one may extend the result of Theorem 11.11 to the parameter regime in which 

113 * 

b = 0{£i). Let b = £i (3 and tj = e* where /3 and 7 are 0(1) with respect to e. Then, the perturbation 
terms in the V 2 component of the non-autonomous system become: 


r 2 



(cos(r|u;f2) cos(r202,o) - l) 


^ sin(r|t<;t2) sin(r202,o) 

^2 


Here, the even terms are 0 {r2) as r 2 —)• 0, so that one may proceed with a similar Melnikov calulation 
as above, and the odd terms again do not contribute to leading order in the Melnikov calculation. 

We further note that a blow-up and Melnikov computation similar to that just presented for the FSN 
I points may be done for the folded nodes and folded saddles, and this gives the location of the maximal 
canards as 

0n,s{k) ~ 2kTr ± cos“^ ^- - - - — h 0{b)^ , /c G Z. 

See also equation (42) in |]2]. 


4 Loci of the Torus Canards and Their Folds for Intermediate- and High- 
Frequency Forcing 


In this section, we study system o in the intermediate-frequency regime with uj = y/eQ and = 
0(1), as well as in the high-frequency regime with oj = 0(1). We prove Theorem |1.2[ demonstrating 
that system ([T]l possesses a family of torus canards in the intermediate-frequency regime, in between 
the two fold curves 0 of these torus canards. The central methods used in the proof are geometric 
desingularization -in which we use a cylindrical blow-up of the fold curve rather than a spherical blow¬ 
up as used in the previous section- and a Melnikov calculation to identify the parameter values for 
which the torus canards exist. After Theorem 1.2 is established, we also prove Corollary |1.3| for the 
high-frequency regime. 

In the intermediate-frequency regime, system ([T]l is equivalent to 


X = y- f{x), 

y'= £{—x + a + bcosO), ( 22 ) 

0 ' = ^/£n. 
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First, we rectify the fold curve so that it coincides with the Q axis, 


u = X — 1, 


v = y + 


2 

3’ 


a = a — 1. 


Also, we recall 


a = ea, 6 = e/3, 


where ci and ^ are 0(1) with respect to e. This transforms ( |2^ to the following system: 

/ 9 f 

U = V — u — -u , 
o 

v' = — £u + (a + P cos 9 ), 

9' = 


Next, we perform the following cylindrical blow-up transformation: 

u = ru, V = e = r^e, 


( 23 ) 


(24) 


which transforms the circle of fold points into a torus. (This contrasts with the spherical blow-up of the 
FSN I point in the previous section.) Append the trivial equation e' = 0 to system ( |2^ and let /r > 0 
be sufficiently small. For each 0 G and for all non-negative values of the system parameters, the 
coordinate change is a map from B := x [—/r, into M^. We examine the vector fields induced by 
( [ 2 ^ in fwo useful coordinafe charfs: fhe enfry-exif chart (or phase-direclional chart) Ki = {v = 1} and 
fhe rescaling chart K 2 = {e = 1 }. 

In chart Ki, fhe coordinates are 


u = riui, V = rf, £ = r\ei. 


Selling F(ui, ri, ei, 0, d, /3) = —ui -h £iri(d + /3 cos(9)), we find fhaf fhe system in chart Ki is 


ui ei 


ui = 1 - 2 

n ei 


(25) 


9 = 

£1 = -sj F, 

where we have desingularised the vector field by rescaling fhe lime variable by a faclor of ri and recycled 
fhe overdol. In fhe phase space of ( |25| ), fhe hyperplanes {ri = 0} and {ei = 0} are invarianf. In addilion, 
for every (a, /3), fhere is an invarianf line 


4 = {(ui,ri,ei) = (rti,0,0)} 

on which fhe dynamics are governed by ui = 1 — uf and 9 = 0. Moreover, for every (ci, (3), fhe poinls 

Pa = (1,0,0) and pr = (-1,0,0) 

are allracling and repelling fixed poinls, respecfively, on iu, and Ihey have Iwo-dimensional center man¬ 
ifolds Na^i and 1 in fhe half-space £i > 0. 

In order lo eslablish fhe exislence of fhe lorus canards, we now hook up fhe dynamics observed in 
charf Ki fo fhose in charl K 2 . In charf K 2 , the coordinates are 

u = r 2 U 2 , V = r|u 2 , e = r^, (26) 

and these coordinates are related to those of chart Ki via the following coordinate transformation: 

-1/2 -1 1/2 

tt2 = til , V 2 =£i , £2= n e/ , 
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where ei > 0. 

In chart K 2 , the system is 


U2 = V2-ul- -r2ul, 


V2 = —U2 + r2 (a + /3 cos(nt2 + ' 


( 27 ) 


where we have also rescaled time by a factor of r 2 in order to desingularise the vector field (with t 2 
denoting this rescaled time variable), recycled the overdot again now to denote the derivative with respect 
to ^ 2 ^ and written the system as a non-autonomous system. For reference, we emphasize the relation 
^/s = r 2 . We now show that system ( |27| ) possesses a special family of homoclinic orbits, connecting the 
point at infinity to itself, which implies that the orbits connect the points pr and pa identified in charf Ki. 


These orbifs correspond fo singular forus canards of fhe original sysfem ( |22l ). 
The unperturbed problem associated fo sysfem ([27l) is given by 


U 2 = V 2 - U 2 , 

V 2 = -U 2 , 

which is fhe same as As shown in fhe previous secfion, fhis unperfurbed sysfem is Hamiltonian 

H{u2, V 2 ) = (^l -V2-^ 

Along fhe level sef T := {H = 0}, which is fhe separafrix befween bounded and unbounded solufions 
(see Figure]^, fhe solufions are given explicifly by 


U 2 ,v{t 2 ) = 


1 

^ 2 ,r(f 2 ) = f- 2 - 


In fhe language of dynamical sysfems, if is a homoclinic orbif fo infinify. 

Wifh fhe above informalion abouf fhe unperfurbed sysfem in hand, we now fum fo show fhaf F 
persisfs for sufficienfly small values of r 2 in ( [27| ). We use a sfraighfforward generalization of Proposition 
3.5 of 136)1, where we nofe fhaf fhe perfurbafion terms there are strictly autonomous, whereas here the 
perturbation terms also include a small-amplitude, time-periodic function, and a compactification of the 
phase space can be used. Moreover, we observe that the parameter there, A 2 , is also treated as being 
a small variable via the linear scaling of A with r 2 , whereas here we have chosen instead to scale the 
parameters a and b with e from the outset and to treat a,/) = 0(1) as parameters. In this manner, r 2 is 
the only small variable in the analysis here. 


The splitting distance between the manifolds Aa ,2 and Ar ,2 for system ([27) is 


D{r2) = dr^r2 H-. 


Here, the dependence of the Melnikov function on the system parameters is implicit. We find 


dr 2 — 


L 


V HI 


— 00 

e 


_1,,3 
3 “2,r 


a + /3 cos(f 2 f 2 + do) 


dt 2 




V 8 


- + a + /3e 2 cos{9o) 


(28) 


where fhe lasf ferm in fhe infegral was evaluated by using 003 ( 2 ;) = Re(e®^), completing fhe square on 
fhe exponenfial, and shiffing fhe confour in fhe complex plane. Hence, reverting fo fhe given parameters, 
we see fhaf fo leading order fhe spliffing disfance is 




a — 1 


+ 



cos(9o) 


(29) 
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Therefore, for each b satisfying the hypotheses of the theorem and for e small enough, the simple 
zeroes of the Melnikov function are given in the (a, fl) plane by 


E 

0 = 1-- — be~^ cos{9o). 
8 


(30) 


This formula, which is exactly Q in the intermediate-frequency regime, gives the parameter values for 


which system pT] ) has a one-parameter (Oq) family of persistent homoclinic orbits, and these persistent 
homoclinic orbits of are the torus canards of ^22\ . 

Also, as a direct corollary, we observe that the envelope of the family of torus canards is given by 


e , / 


(31) 


which is precisely formula Q. This completes the proof of Theorem |1.2| 

The proof follows by extending, in a straight 


To conclude this section, we prove Corollary 1.3 


forward manner, the above analysis of the persistent homoclinics and the folds of the torus canards in 
the proof of Theorem |1.2| for the intermediate-frequency forcing regime to the high-frequency forc¬ 
ing regime. In particular, in the high-frequency regime, oj = 0(1), which corresponds to taking 
n = 0{l/^/e) in the above analysis. Following exactly along the above calculations, we see that 
the geometric desingularization method yields the same equations ( [27l ) in chart K 2 , but now the small- 
amplitude time-periodic forcing term has high-frequency 12 = 0{l/^/e). Hence, the suitable version 
of the Melnikov theory is that for rapidly forced systems, and the splitting distance along T is again 
given by ( |3T] ), which is now exponentially small in e, since 12 = 0{l/^/e). See for example lfT2ll . This 
completes the proof of Corollary |1.3| 

The result of this Corollary for the high-frequency regime also agrees well with the results obtained 
from numerical simulations. In Figure]^ for oj = 0(1), we present a computation of the distance 
between the two folds of maximal primary canards as a function of e. We gathered the control points 
obtained for various computations for eleven fixed values of e, decreasing from 3 • 10“^ down to 8 • 10“^ 
and plotted them on a logarithmic scale. The hyperbolic shape of the resulting curve confrms that this 
distance is exponentially small in e as e tends to 0 . 



Figure 7: Width of the canard region in the uj = 0(1) regime, as a function of e for & = 0.01. 


5 Secondary Canards 

Having established the existence of the primary strong canards and their folds, we now turn our attention 
to the secondary canards of the folded nodes of ([T]l, which exist in the low-frequency forcing regime oj = 
e io. By definition, secondary canards lie in the transverse intersections of the invariant slow manifolds 
and S^. A representative example of these manifolds and their intersections (i.e. the secondary 
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canards) is shown in Figure These manifolds are computed from curves of initial conditions traced on 
the attracting and repelling sheets, respectively, of the critical manifold S, up to a cross-section at fixed 
angle 9 corresponding to the maximal torus canard lfT3l [T4ll : 


^77, — ^ 9ri. — COS 


_i /1 — a — e /8 


+ 0{b) 


In Section 5.1 we study the folds of the secondary canards and investigate how they change in the (ce, a) 



Figure 8: Attracting (red) and repelling (blue) slow manifolds of system Q for a = 0.9935, b = 0.01, w = 0.005, 
and e = 0.05, together with a stable MMO with 1 LAO and 6 SAOs found for the same parameter values by direct 
simulation. Bottom row; invariant slow manifolds (red) and (blue) of system 0 in the cross-section 
The secondary canards are identified as the intersections (black dots) of S'® and S®. 


plane under variation of the forcing amplitude b (analogous to the folds of the primary canards). We then 
investigate in Section [S!2| how large-amplitude oscillations can grow from small-amplitude oscillations. 


5.1 Continuation of Secondary Canards 


As shown in Figurej^ the invariant slow manifolds S® H and S® H spiral around one another, which 
is typical of a folded node. More precisely, let fj, := X^/Xg, |A 7 „| < |As|, denote the eigenvalue ratio of 
a folded node, regarded as an equilibrium of the desingularised system ( [T^ . Provided e is sufficiently 
small and /r is bounded away from zero, the total number of (primary and secondary) maximal canards 
is Smax + 1, where 


'^max 


L—J- 

^ 2/i 


and [-J denotes the floor function. In particular, a persisfenf branch of secondary canards bifurcales from 
fhe weak canard in a franscrifical bifurcalion for odd integer values of |[55l . 
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Remark 3. The A:* secondary canard exhibits k small oscillations about the weak canard for k = 
1,2,, Smax — 1- These small oscillations are localized to an 0(^/e) neighbourhood of the folded 
node Il55ll5^ . Moreover, trajectories on situated between jk-i and 'y/-, k = 1, 2,..., Smax execute 
k small oscillations about the weak canard, where 70 and jsmBLK correspond to the primary strong canard 
and primary weak canard, respectively. 

By tracking the resonances = 2A: + 1, fc = 0,1, 2,..., we can follow (in the singular limit) the 
locations in the (uJ, a) plane where the secondary canards are born. Figure|^shows an example for 6=1. 
The non-singular (co, a) plane shows that only the folds of canards corresponding to the FSN I and the 
degenerate folded node extend into the intermeditate frequency regime. All other branches of folds of 
canards are restricted to the low-frequency regime w = 0 (e). 

Remark 4. Note that the resonance curves in Figure bear no resemblance to the curves of folds of 
secondary canards in Figure]^ This is to be expected since 6 = 0(^/e) in Figure]^ which implies that 
/r = 0 (\/e) and the folded node theory does not apply. 




Figure 9: Resonance curves for 6=1, and (a) £ = 0 and (b) e = 0.01. In (b), the theoretically computed 
curve of maximal canards for the FSN 10. are shown in red. Inside this envelope, there are two black resonance 
curves. Both correspond to the numerical continuation of folds of maximal canards. The outermost black curve 
corresponds to the FSN I (i.e. where /r = 0). Note that 0 breaks down when oj is no longer 0(e). The inner black 
curve corresponds to the maximal canard of the degenerate folded node (i.e. where /j, = 1). The inset shows the 
numerical continuation of the folds of canards corresponding to /r = 1 and /r = 1/3. In particular, for the fJ- = 1 
resonance, we compare the numerically computed result (blue) and the theoretical result obtained in ( [3^ (red), 
and find that there is excellent agreement away from the FSN I boundaries. 

For the degenerate node (/i = 1), a Melnikov computation similar to that in Sections]^ andshows 
that the locus of the primary maximal canard of the degenerate folded node in the (uJ, a) plane is 

a = l-|±yi,= ^exp(-ljiS^), (32) 

which holds provided uJ = 0(1), = 0{e) and y^6^ ~ (1 ~ = 0(i/e). The inset 

of Figure shows that there is excellent agreement between this theoretically computed curve and the 
curve obtained from numerical continuation. The deviation between theoretical and numerical results 
for this degenerate folded node maximal canard starts to become significant when the degenerate node 
branches approach the FSN I branches. We note an important implication: all secondary canards due to 
folded nodes are restricted to the region of the {uj, a) plane bounded by w = 0, the locus of the folds of 
maximal canards of the FSN I and the locus of the maximal canards of the degenerate folded node. 

As was the case for the numerical continuation of the folds of the primary canards and the folds 
of the torus canards, the numerical continuation of the maximal secondary canards is done by solving 
families of boundary-value problems and computing branches of such solutions using pseudo-arclength 
continuation. Along these branches, a number of fold points can be detected, and then the curves of folds 
of secondary canards can be continued in two parameters. For a representative set of parameter values. 
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the folds of the first, second, ..., tenth secondary canards (i.e., with respectively one, two, ..., ten loops) 
are shown in the (cj, a) plane in Figure]^ The outermost envelope in Figure |^are the curves of folds of 
primary canards. The rightmost path enclosed by the fold curve of the primary strong canards represents 
the fold curve of the first (1-loop) secondary canard, and each successive curve to the left represents a 
family of folds of secondary canards with one additional loop. These branches of folds of secondary 
canards emanate from the FSN I points at a = 1 — | zh 6 . As a; is increased, the corresponding pairs of 
n-loop branches come together at turning points. 

It is also useful to examine projections of the secondary canards onto the (x, y) plane. In Figure [l^ 
we show the first three maximal secondary canards, with respectively, one loop (yellow), two loops 
(red), and three loops (blue). The highest loops of the 2-loop and 3-loop maximal secondary canards are 
observed to lie extremely close to the single loop of the 1-loop maximal secondary canard. The same 
holds for all of the higher-loop secondary canards, as well. Also, the second loop of the 2-loop canard lies 
inside the first loop, and it lies extremely close to the second loop of the 3-loop canard. In addition, for 
even smaller values of the forcing frequency uj, the y-intercepts of the return jumps increase. Moreover, 
these ^-intercepts diverge to oo in the limit a; —)■ 0. In fact, in this limit, the maximal secondary canards 
collapse onto the primary strong canard, consistent with the observation that the branches of maximal 
secondary canards emanate from the same FSN I points as the primary canards do. 



Figure 10: Projection of the 1-loop (yellow), 2-loop (red), and 3-loop (blue) secondary canards of 0 shown 
for UJ = 0.01, £ = 0.01, b = 0.01, and a = 1.0034909031 (1-loop), a = 1.0034909029 (2-loop), and a = 
1.0034909029 (3-loop) shown in the (x, yj-plane. 

We now investigate what happens to the invariant slow manifolds near the turning points (recall 
Figure]^ of the fold curves of secondary canards. In Figurewe show one such fold curve and take 
four values of uj, for a fixed value of a, near the turning point which marks the largest cj-value of this 
curve (top panel). For each value of uj, we compute and S!^ up to a fixed cross-section, following the 
procedure described above. Then, the intersection curves of both manifolds in the fixed cross-section are 
shown for each w-value in the four bottom panels of Figure Each time the fold curve of maximal 
canard solutions is crossed, two intersections of the attracting and repelling slow manifolds disappear or 
are created. This is illustrated in each of the transitions shown in panels (l)-(4). 

5.2 Growth of LAOs from SAOs in the Secondary Canards 

Along the continuation of the secondary canards, an orbit segment can ‘grow’ an LAO. This occurs in 
regions where the repelling slow manifold spirals backwards instead continuing to spiral inwards towards 
the weak canard. In Figure [T^ we show the repelling slow manifold for uJ = 0.3 in the cross-section 
Tin', the direction of spiralling changes three times along this portion of the slow manifold, at the points 
labeled (hi), (b2), and (b3) in frame (a). Each direction reversal corresponds to the orbit growing a 
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Figure 11: Evolution of the intersection points (black dots) between the curves representing and in a fixed 
cross-section through the folded node. Here, a = 0.99641, b = 0.01, and e = 0.05. The values of w are (1) 0.049, 
(2) 0.0505, (3) 0.0517, (4) 0.0525, as also labeled on the horizontal axis in the top frame. In the transition between 
frames (1) and (2), the lower two intersection points disappear. Between frames (2) and (3), two intersection points 
are created. Then, two disappear in the final transition shown, from (3) to (4) 


large-amplitude oscillation. In panels (bl), (b2), and (b3), we show the profiles of the computed solution 
segments at these events. Note that the first fold encountered (starting from the center of the spiral and 
going outwards) corresponds to the orbit segment having one SAO grow up to the size of an LAO (this 
occurs at the second fold), while the other SAOs stop growing in size. 

Remark 5. An important consequence of studying the curves of maximal canards and maximal torus 
canards in the parameter space of ([T]l is that they serve as the boundaries between different dynamic 
regimes of O- As highlighted briefly by fhe graphical summary in fhe {oj, a) plane shown in Figure|^ fhe 
fvdP equafion ([T]l exhibifs small-amplifude oscillafions (SAOs), large-amplifude or relaxafion oscillafions 
(LAOs), and mixed-mode oscillafions (MMOs). The SAOs are fhe --periodic solufions generafed when 
an affracfing equilibrium of fhe unforced vdP equafion (i.e., 6 = 0) is subjecfed fo a small-amplifude 


21 























Figure 12: Repelling slow manifold reversing the direction of spiralling as successive LAOs appear in the orbit. 
Panels (hi), (b2), and (b3) show the solution profiles corresponding to three such events. Here, a = 0.9935, b = 
0.01, uj = 0.3, and e = 0.05. 


periodic forcing of frequency to (Figure [T3]^a)). The LAOs occur when the equilibrium of the planar 
vdP equation sits on the middle branch of the cubic-shaped nullcline and the attractor of the system 
is a relaxation oscillation that alternates the trajectory between the outer branches of the cubic (Figure 
[T3]^f)). The MMOs feature SAOs superimposed on large-amplitude relaxation-type oscillations. Figure]^ 
shows that the MMOs become more robust for low-frequency forcing, just as was observed in numerical 
simulations of a rotated van der Pol-type model in f2l ■ 

6 The fvdP is a Normal Form for a Class of Systems Near Torus Canard 
Explosions 

In this section, we prove that, under a number of natural conditions, a slow/fast system with two fast 
variables and one slow variable, which is subject to time-periodic forcing and for which the fast system 
possesses a generic fold of limit cycles, is equivalent to a system in which the fast component is given to 
lowest order by the forced van der Pol system ([T]). We consider systems with two fast variables and one 
slow variable of the following form: 

X = fi{x,y,z) 

y = f 2 {x,y,z) (33) 

z = eg{x,y,z), x,y,z£R. 


The fast system is 


X = fi{x,y,z) 
y = f 2 ix,y,z) 


(34) 
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Figure 13: Transition from SAOs to MMOs to LAOs in the forced vdP equation ([T]). The attractors of 
([T]) are shown for e = 0.01, u = 0.01, and (a) a = 1.01, (b) a = 0.9999, (c) a = 0.999, (d) a = 0.994, 
(e) a = 0.991, and (f) a = 0.987. 


in which z is a parameter, and we make the following hypotheses about the fast system: 

(HI) For z = 0, there exists a non-degenerate periodic solution F; and, here by non-degenerate, we mean 
a periodic solution with finite period. 

(H2) The Floquet multiplier of this periodic orbit at z = 0 is one. 

We prove the following theorem: 


Theorem 6.1. Given (HI), (H2), and a non-degeneracy assumption (see below), system ( |33l ) is 
locally (in a small neighborhood of F) orbitally equivalent to 


p = z- + 0{p^) + 0{e) 

e = 1 (35) 

z = £g{p,6,z), p,z,9£R, 


where g is 27r-periodic in 9, and where z and p are small. 

The fast subsystem of (^i has very similar dynamics to that of ([T]l, however the slow system of the 
general systems may be much richer than the slow system of ([T]l. We make a more detailed comparison 
between the full systems later in this section. We first analyze the fast systems. 

We begin the proof of Theorem |6.1| with some preliminary transformations. On the basis of (HI), one 
may rectify the flow of ( |34l ) so that the periodic orbit becomes the unit circle. Next, using the coordinates 

X = (1-|-r) cos(0) and y = (1-|-r) sin(0), 


one may transform (pd]) to 


r = Mr,9,z) 
0 = f2{r,9,z), 


(36) 
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where /* is 27r periodic in 9, /i(0, 9, 0) = 0, and /2 7 ^ 0 in a neighborhood of (0, 9, 0). Also, one may 
scale the time variable so that (|3^ becomes 


f = F{r, 9, z) 
9 = 1, 


(37) 


with F = /i// 2 - This is a useful formulation of the fast system, and we work directly with this through¬ 
out the proof. 

In order to analyse the dynamics of this system for small values of r, we expand: 


F{r, 9, z) = 'iIjo{9, z) + rV'i(0, z) + r‘^ip2{0, z) + ... + z) + 

Condition (H2) is now equivalent to 


p27r 

I ' 


V’i(0, 0)d9 = 0, 


(38) 


(39) 


because 


where the average of F is defined as 


dF 

^(0,0) = 0, 


1 


r-27r 


F'{r,z) = — F{r,9,z)d9. 

Jo 

Also, to analyze the dynamics for small values of r, it is useful to introduce a new variable p by 


r = (/>o + /oe‘^1 + /9^(/>2 + ... + p^'^4 >n, N >2. 

where (j)i are functions of {9, z) to be determined. 

The following lemma lies at the heart of the proof of this theorem: 

Lemma 6.2. There exists a choice of functions (pi, i = 1,..., A^, in the coordinate change (|4^ and 




(40) 


function G{p, z) such that, in the coordinates {p, 9), the fast system ( [37] ) has the form 

p = G{p,z) + 0{p^"^^) 

9 = 1. 

Remark 6. Hypothesis (H2) is not needed for this lemma. 


(41) 


Proof. The proof is split into three steps. First, we consider the simplest case in which N = 1 in (40). 


Next, we prove the Theorem for N = 2, and finally we prove if for general in (401. 


Step 1. Wifh A^ = 1, fhe relevanf fransformafion ( |40| ) befween r and p is: 

r = (j)Q + . 


(42) 


Differenlialing ( |42| ) wifh respecf fo time, subsfifufing in ( |T7] ), and Taylor expanding F(i^o + , 9, z) 

abouf pQ, we obfain 


P = {-(po,e + F{(t)Q, 9, z))e + pi-pi^ + T;(</>o, G, z)) + O(p^). 


(43) 


Hence, based on fhe form of fhe fwo ferms in parenfheses in fhe righf member of fhis equation, we are 
nafurally led fo sfudy fhe following system of fwo ordinary differenfial equafions wifh fwo unknown 
parameters: 


cPi,^e = F{Po,9,z) + \oe^^ 
(t>i,9 = Fr{(j)0,9,z) + Ai 


(44) 
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In this manner, proving the result for = 1 is now equivalent to finding Aq and Ai for which there exist 
(j)Q and (f)i which are 27r periodic and which satisfy the equations ( |4^ . 

According to general theory of differential equations, for every pair (Aq, Ai), there exists a unique 
solution ((/iQ, (/>i) of ( |4^ satisfying the initial conditions </>o(0) = 0 and (/>i(0) = 0. Such solutions must 
satisfy the integral equations. 


p9 p6 

(paid) = / F{(po{u),iy,z)diy + Xo / 
Jo Jo 

= f + XiO. 

Jo 


(45) 


Also, in terms of these integral equations, the conditions for periodicity are 


p27T ^27r 

0 = / F{(j)Q{u), v, z)di'+ Xq / 

Jo Jo 

p27V 

0 = / Fr{4>oii^),’^, z)di'+ 27rXi. 

Jo 


(46) 


Based on the above formulation of the integral equations, it is useful to define : M x 
as follows: 

, \ t ?fo( 2 ,Ao,Ai) F{4>o{u),iy, z)du + Xq ^ 

^ V -^o> ^i) / A Jo"" Fr{(poi’^),i^,z)du+ 2 ttXi j 


(47) 


Now, by the assumptions (HI) and (H2) (see (|3^), Aq = 0, Ai = 0, and z = 0 is a solution of (451 
with (/)o = 0 and i^i o = Jq Fr{0, u, 0)du. We now verify that the assumptions of the Implicit Function 
Theorem are satisfied to show that there is a branch of nontrivial solutions emanating from this trivial 
solution. In particular, we verify that D'H{0,0,0) is non-singular, by showing that detDTf (0,0,0) = 
( 2 vr) 2 . 

From the definition, we see that 77 i,Ai (0,0,0) = 27r. We will now show that 77 o,Ai (0,0,0) = 0 and 
77o,Ao (0) 0) 0) = 27r. To show that 77o,Ai (0,0,0) = 0, we start with 


r2n 


77o,Ai ( 0 , 0 , 0 ) = / FriO,e,O)cPo,xMd0- 


Then, using ( [44l ), we obtain 

'Ho,Ai(0,0,0) = ■^(p0M{^)d0 = (^o,Ai(27r) - ^o,Ai(0). 

By 3,ssumptiori (3.bovc forinuls. ( |45| ))) (^q(O) ^ 0 3 .n(d (^i(O) ^ 0 for Aq 3 .ncl A)[. Hcncc, <?^o,aAO) = 0 
for j = 1,2. Therefore, we conclude from ( [44l ) that ((>o,Ai = 0, so that this entry of the Jacobian of 77 
vanishes, as claimed. 

Next, we show that 77o,Ao(0, 0,0) = 27r. By an argument similar to the above. 


77o,Ao( 0,0,0) = ())o,Ao(27r) - (/>o,Ao(0) = ())o,Ao(2vr). 


Therefore, from (44 1 , we obtain 


A 

dO 


i4>o,\o) = Fr{0,6,0)4io,Xo +e'^T 


(48) 


Observe that, for Aq 


0 and Ai = 0, we have dcpi/dO 


Fr{0,9, 0). Hence, (481 is equivalent to 




d4>i 

de 


4’o,\o + , 


(49) 
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from which it follows that 


Finally, since cpo^Xo (0) 


A. 

dB 



= 1 . 


0, it follows from (^i that = 271. Therefore, 


det(Z)'H(0,0)) = (27r)2, 


and by the Implicit Function Theorem there is a branch of periodic solutions (j)o and cpi for each (Aq, Ai) 
sufficiently small, emanating from the trivial solution. This completes the proof of the Lemma for the 
case = 1. 


Step 2. We now show that the lemma holds for = 2 in ( |4()| ). All quantities are expanded up to and 
including p^. For the vector field F, we have 


F(r,e,z) = F{4>o,9,z) + pFr{(l)o,6, z)e'^^ + p‘^{Fr{4>o,d, z)(l )2 + ^Frri(po,d, 

Also, differentiating ( |40| for A^ = 2 with respect to t, we find 

f = + 2pcj)2) + (po^d + + P^(p2,e- 

Hence, combining ( [SO] ) and ( |5T] ), we get 

p{l + 2 p(j) 2 e~'^^) ={-(po,e + F{4’o, 0, z))e~‘^^ + pi-cpi^e + Fricpo, 6, z)) 

+ p'^{-<t) 2 ,e + Fr{4>o-,0,z)<t)2 + ]^Frr{<t)o,9, . 


(50) 


(51) 


(52) 


Now, after multiplying both sides of ( |52l ) by (1 + 2pcj)2e~'^^)~^, we examine the structure of the 
terms at each order of p^, p^, and p^. This suggests that we analyze the following system of differential 
equations: 


<Pofi = F{(Po,e,z) + \oe^^ 

(pifi = Fr{(j)o, 0, z) - 2Ao02e“'^^ + Ai 

(p 2 ,e = Fr{(t)o,0,z)(j)2 + ^Frri(l)o,9,z)e^'^^ - 2Xi(j)2 + A 2 e'^L 


(53) 


The first equation here is equivalent to the first equation in ( |4^ ; the second has an additional term due 
to (1)2', and, the third is new. If we can find a branch of nontrivial solutions of this system, then we can 
transform the general fast system into the desired form up to and including 0{p^). 

As above, we look for 27r periodic solutions. However, before extending the definition of l-L, we 
rewrite the third component of ([5^ as follows: 


(54) 


A = \Frr{(Po, 9, z)e^^ + 2Xo4e-^^^ - 3Xi(j)2e-^^ + Aa- 


Now, we define V. : 


in a manner similar to that employed in Step 1: 


F-oiz, Ao, Ai, A2) 

Tf(z, Ao, Ai, A2) = ( T-li{z, Xq, Xi, X2) I = 

F,2{z, Ao, Ai, A 2 ) 

^ Jq"" F{(j)o{v), V, z)di' + Ao Jq"" 

z) - 2Ao</>2(z^)e“‘^i(''))di^ + 27 rAi 

V Io"'ilFrr{(j)o{D), D, + 2Ao(()2- 3Xi(j)2{y)e~‘^'^^^^ + X2)dv 


(55) 
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We analyse % in much the same manner as in Step 1. Observe that, at z = 0, we have Aq = 0 and 
Ai = 0. However, A 2 is in general given by 


A 2 — 





and not zero. 

We now show that the off-diagonal elements of the Jacobian of Ti vanish. First, an argument similar 
to that used in Step 1 shows that (0; 0) 0; -^2,0) = 0- We also need to show that Hq ,A2(0, 0,0, A2,o) = 
0 and Hi,a 2 (0; 0) 0; ■^2,0) = 0. To this end, we prove that ^o,A 2 (0) 0) 0) "^^2,0) = 0 and 4’i^x2 (0; 0, 0, A2,o) = 
0, because then the identities ?fo,A2 (Oj Oj 0) -^2,0) = 0 and i,a 2 (Oj Oj 0) "^2,0) = 0 follow in a straight¬ 
forward way. We carry out the proof for 0o,A2i the argument for 0 i,a 2 similar. Differentiating the first 


equation in (531 with respect to A 2 and using the fact that Aq | 2 =o = 0, we obtain 


d4’o,X2\z=o 


= F;(0,6I,0)(/)o,A2|2=o- 


(56) 


The claim now follows from the assumption (/)o(0, z, Aq, Ai, A 2 ) = 0. 
Based on the above analysis, it follows that 


det{Dn{0, 0, 0, 0)) = Ho,Ao(0, 0, 0, 0)'Hi,Ai (0,0,0,0)^2,A 2 (0, 0, 0, 0), 

and an argument similar to the one used in Step 1 shows that this determinant is nonzero. In par¬ 
ticular, ?fo,Ao (05 0) 0) 0) Tfi^Ai (0,0, 0, 0) are both 2tt by a similar calculation. To show that also 
^ 2 ,A 2 (0) 6) 6) 6) = ^TT, we differentiate the third component in ( |55l ) to obtain 

r2n 

H 2 ,A 2 ( 0 , 0 , 0 , 0 ) = / dv = 2'n. 

Jo 

Hence, we may again use the Implicit Function Theorem to conclude that there exists a branch of non¬ 
trivial solutions of ( [5^ , and the system may be put in the desired form up to and including terms of p^. 
This completes the proof of Lemma [6^ for the case N = 2. 


Step 3. In this third and final sfep of the proof, we show that the lemma holds for general N in ( |4^ . We 
begin by writing 


N 


F{r,e,z) = ^p>Fj{(j)o,(f)i,... ,4>n,0,z) + ( 9 ( p ^+^), 
i=o 


(57) 


where N > 2 is a natural number and r and p are related by formula ( |40l ) associated to this choice of 
N. For each j, the functions Fj are complicated expressions involving (j)Q, (f)\, ..., cjij. To simplify the 
notation, we write = [cpo, ^j). We will give a more precise description of the functions Fj 

below. The equivalent of (|5^ is now 


N 


p(l + ^(p^ ^4>ie ={-(t)o,e + F{(t)Q,9,z))e + p{-(t)ifi + Fr{(j)Q,9,z)) 


1=2 


N 


( 58 ) 


j=2 


Let ao = 1 and for each j > 1 set 


aj{4>0 ,... ,9)n) = 


j\ ^Vl + Er=2^p' 
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Note that, for given j, aj depends on (j)o, (j)i, ..., (j)j+i. We also write 

Eo = {-(l)o,e + F{cl)o,e,z))e-^^ 

El = -4>i^e + Er{(po,0,z) 

E2 = {-cj)2,e + F2{^2,0,z))e-<^^ 


( 59 ) 


En = i-<pN,e + Fn{<^n, 9, z))e-<^K 


It follows that (p8|) maay be written in the following compact and insightful manner: 


N / I 

z=o \i=o 

We now define the set of equations 

Eq = Ao 
El + uiEq = Ai 
E 2 + cci-Ei + a2Eo = A2 



N 

'^ajE]\f_j = Aat. 
j=0 

This enables us to rewrite l59l as follows: 

Eq = Ao 

El = Ai — ctiAo 

E2 = X2 — aiAi + (af — q; 2 )Ao 

N 

En = Pj^N-j, 

j=0 


(60) 


(61) 


(62) 


where /3o = 1 and f^i,... ,(5 n are coefficients depending on oq, • • ■, oat. Moreover, fOj depends on 
aO) • • •) ctj only. Therefore, we have arrived at the following system of differential equations: 


ct>o,e = F(cl>ii,e,z) + \oe^^ 

4>i,e = Fr((po, 9, z) + Ai - aiAo 

<p 2 ,e = F2(^2',9., z) + (A 2 — aiAi + (af — a 2 )\o)e'^^ 

(63) 


N 


4>N,e = Fn(^n,9,z) + j '^l3j\N-j j 

vi=o 


which is the analog for general N of the systems of differential equations ( [4^ for A/" = 1 and ( [53] ) for 
N = 2. 


Before defining "H, we rewrite ( |63) in a manner similar to that which was used above to rewrite ( |53| ) 
(recall also (|5^). Noting that 


Fj{^j,9, z) = Fr{(l)o, 9, z)(j)j + 7^(4>J_l), 
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we replace the jth equation in ( | 6 ^ by 






) = - (Ai - • 


(64) 


We will now define the function in a manner analogous to ( |47l ) and (^i used in STeps 1 and 2, 
respectively, for the cases N = 1 and N = 2. We let 

((/>o(6*, Z, Ai, . . . , Aat), 0i(0, Z, Ai, . . . , Aat), . . . , (/)Ar(0, Z, Ai, . . . , \ n )) 


be the solutions of ( [6^ depending on the parameters z and Aq, • • •, A^v that satisfy the initial conditions 
(?!)j(0, z, Ai,..., Aat) = 0, y = 0,1,..., A^. Further, we let 1-Lq, T-Li, and 4^2 be defined as in Step 2, and 
let 

Tijiz, Ao, Ai,..., Aat) = J — (Ai — aiXo)4>je~'^^ + ^ du, y = 2,3,. 

(65) 

We first argue that we can solve the set of equations 'Hj{0, Aq, • • ■, Aat), j = 0,1,..., N, for a unique 
A^-tuple (Ao, 0 ) Ai_o, • • •, AAr,o)- Note that Ao,o = 0 and Ai_o = 0, by the same analysis used in Step 2. 
Hence, 

j i-2 

^ ^ /5fcAj—fc|z=o — ^ fdkXj—k\z=o- 
k=0 k=0 

We argue by induction. Suppose that Ao,o, • • •, Aj_i^o and the corresponding 4)0,..., 4>j-i are deter¬ 
mined. Note that Ayo must satisfy 


r- 27 r 


i -2 


Ajp - - 


7^(<h,■_l)e-'^l + J]AA 


j-l,0 


dv. 


( 66 ) 


1=1 


Since the right member of ( [661 ) depends only on Aq, • • •, Aj_i and ^o,..., (pj-i the value of Xj^ is 
uniquely determined. Similarly, knowing Aj_o^ we can solve ( |64| ) for 4>j using ( |64| ) due to the fact that 
Aq = Ai = 0, so that right member of ( |64| ) depends only on Aq, • • •, Aj_i and 4>o, ■ ■ ■, (4j-i- Hence, 4>j 
is uniquely determined. 

We now prove that Ao,o, • ■ •, Aat.o) N non-singular. First, we prove that 4>jM U=o = 0 

for any j G {0,..., — 1} and k £ {j + 1, ■ ■ ■, N}. The argument for j = 0 and j = 1 is analogous 

as in the case of n = 2. For general j, we proceed by induction, assuming that the claim holds for 
0,1,..., y — 1. The argument is, again, similar to that used in the case N = 2. We differentiate (Ml with 
respect to Xk and use the induction assumption, the fact that is independent of 4>i for any l> j + l, 
and the fact that Ao,o = Ai^o = 0 . This gives 


Ao,o,..., Aat.o) = 0. 


(67) 


By assumption, 4>j{0, z, Xq, ■ ■ ■, X^) = 0. Hence, the claim follows. Now, differentiating ( [65] ) and 
using a similar procedure, we obtain T-Lj^x^i^, Ao, 0 ) • • ■ j AAf,o) = 0 for y G {0,..., At — 1} and k G 
{j + l,...,At}. 

It remains to prove that Hj^Xj (0, Ao,Oj • • •, Aat.o) 7 ^ 0 for j G {0,..., At}. The proof for j = 0 and 1 
is as above. Let j > 1. Again by differentiating ( |^ , now with respect to Xj, and arguing analogously 
as above, we obtain 


r2TT 


77j,Aj (0, Ao,o, • • •, AAr,o) = / dv = 2'K. 

Jo 


( 68 ) 


We can now apply the Implicit Function Theorem to obtain the functions Xo{z),, Xj(z) as required. 
This completes the third (and final) step of the proof of the lemma. 


..N. 
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QED 


Proof of Theorem 6.1 If we apply the sequence of transformations leading to ( |4T] ) to system ( |3^ , we 
obtain a system of the form 


with 


P — G{p, z) + + 0 {e) 

9 = 1 

i = eg{p,d,z,£). 


N-l 

G{p,z) = 

j =0 


(69) 


and the coefficient functions are as introduced in the proof of Lemma 6.2 In particular hypotheses (HI) 
and (H2) imply Ao(0) = 0 and Ai(0) = 0. We now formulate the additional degeneracy assumptions in 
terms of the functions Xj. We will assume that N > 3. 


A2(0) / 0. 


(70) 


If ( |70| ) holds in addition to (HI) and (H2) we replace the variable z by z = Ao(^;) and perform scalings 
and translations to arrive at ([35]), taking e sufficiently small as necessary. QED 


We now make some remarks relating the normal form equation ( [35] ) derived in this section and the 
forced van der Pol Q- Clearly, the fast subsystems are similar. If in addition to ( |70| ) we assume that 
-^ 3 ( 0 ) 7 ^ 0 and make some assumptions about the signs of the coefficients, then the two fast systems are 
the same to lowest order (up to a simple transformation). The situation with the slow equation is more 
complicated. In systems with time-periodic forcing, we have shown that ([T]) is a normal form. Moreover, 
because forced systems are not generic in the larger class of general slow/fast dynamical systems, we 
expect the dynamics in this larger class to be even richer. In addition, if we assume that p(0, 9o, 0, 0) = 0 
for some 60 G [0, 27r) then we may obtain some canard dynamics. In general these slow/fast systems 
will have folded singularities and associated canard dynamics. 


7 Conclusions and Discussion 

In this article, we have established the existence of a number of different types of canard solutions of the 
forced van der Pol equation ([T]) across the entire range of forcing frequencies a; > 0. Most interestingly, 
we have found numerically that the families of primary maximal canards and maximal torus canards are 
organised along single branches in parameter space. In the low-frequency regime {to = 0(e)), Theorem 
1 1.1 1 demonstrates the existence of the primary maximal canards of the ESN I points and establishes that 
formula ([^ gives the loci of the folds of the primary maximal canards in the (a, b, uj) parameter space 
with b = 0{^/£). In the intermediate-frequency {uj = 0{y/£)) and high-frequency regime (w = 0(1)), 
Theorem 1 1.2 1 and Corollary ] 1.3] respectively, establish the existence of maximal torus canards as well as 
the formulas (]^ and (0, which explicitly give the locations of the folds of the maximal torus canards in 
the (a, 6, cu) parameter space with 6 = 0(e). These maximal torus canards lie precisely in the intersection 
of the persistent critical manifolds of attracting limit cycles and of repelling limit cycles. They are the 
analogs in one-higher-dimension of the maximal headless ducks of the unforced van der Pol equation, 
see for example E FI^ITTI . Moreover, they are similar to the folds of maximal torus canards observed 
earlier in a rotated system of ven der Pol type, see Eigure 5 in ]S]]. 

It was also shown that these analytical results are all representations of the same formula (]^ that 
holds across the entire range of forcing frequencies a; > 0 for the appropriate values of b, and that these 
formulas agree well with the results obtained from numerical continuations over the parameter regions 
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in which they apply. Moreover, in the limit w —)• oo, the torus canards appear to be rotated copies of the 
limit cycle canards that exist in the planar unforced van der Pol equation, and the interval of a values for 
which the maximal torus canards exist shrinks to the value ac(e) = 1 — | at which the maximal headless 
duck solution exists in the unforced equation, recall 

It is worth noting that the analytical results presented here for the torus canards of ([T]l agree with and 
expand upon by the general topological analysis presented in Section 6 of ||3. There, fast-slow systems 
with two fast variables and one slow variable were studied in which there is a torus canard explosion in 
the transition regime from stable periodic spiking (tonic spiking) to bursting, and examples were given, 
including of the Hindmarsh-Rose equation, the Morris-Lecar-Terman model, and the Wilson-Cowan- 
Izhikevich system. In particular, it was shown using topological arguments that there must be a sequence 
of torus canards in these transition regimes in order to satisfy the property of continuous dependence of 
solutions on parameters. The topological analysis presented in ||8l is the analog in one higher dimension 
of the topological analysis first used in fJl to establish the existence of an explosion of limit cycle 
canards in the transition between asymptotically stable solutions and full-blown relaxation oscillations 
in the unforced, planar van der Pol equation. 

In this article, we also studied the branches of the secondary maximal canards, which exist in the 
low-frequency regime in ([T]l. Secondary canards lie close to the primary strong canard for most of their 
lengths, and in addition they make finitely many loops near the bottom of T, recall Figure 7. They are 
indexed by the number of loops and by the height in the y variable of the jumps from the repelling slow 
manifold back to the attracting slow manifold. We showed how the dynamics of these secondary canards 
changes as the parameters change along the fold curves, and we identified the mechanism by which these 
branches turn around well before they get into the high-frequency regime. In particular, the turning points 
correspond precisely to the parameter values at which the fold curve of maximal canards is crossed and 
two intersection points of the attracting and repelling slow manifolds are created (or annihilated). In 
addition, we identified how new LAO segmenfs are added fo fhe secondary canard solufions af poinfs af 
which fhe direcfion of spiralling of fhe repelling slow manifold is reversed, recall Figure 9. 

Finally, we proved fhaf the fvdP equation ([T]l is a normal form for a class of slow/fast systems with 
two fast variables and one slow variable, which possess a non-degenerate fold of limit cycles in the 
fast system and which exhibit the torus canard explosion phenomenon. Thus, the methods and results 
obtained here for Q extend naturally to a large class of slow/fast systems with single-frequency time- 
periodic forcing. 

To conclude this article, we discuss a number of topics related to the canard solutions of the fvdP 
O- First, the fold curves of the primary maximal canards in the low-frequency regime and the fold 
curves of the maximal torus canards in the intermediate- and high-frequency regimes together serve as 
the boundary of the MMO regime in ([T]l, recall Figures and [T3] 

Second, there are many different branches of resonance curves, curves of torus bifurcations, saddle- 
nodes of periodic orbits, period-doubling curves of periodic orbits, and so forth, all of which lie in the 
interior of the MMO region in parameter space. These bifurcation curves play important roles as the 
boundaries between orbit segments with different numbers of SAOs and LAOs. These are the subject of 
ongoing research. 
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A Proof of the Existence of a Torus Bifurcation 


In this appendix, we consider the forced vdP oscillator ([T]) subject to high-frequency forcing (w = 0{1)). 
In this case, y is the only slow variable, and there is no critical manifold. As such, the fast dynamics 
are dominant throughout the phase space, and system ([T]l can be interpreted as a regularly perturbed, 
non-autonomous problem. 

We establish the existence of a torus bifurcation using second-order averaging and equivariant normal 
form theory in Section A.l[ and we present the asymptotics of the torus bifurcation parameter value 
aTB(^, s) in Section A.2 


A.l Second-Order Averaging and a Torus Bifurcation in the Regime cu = 0{1) 

In this section, we perform a standard second-order averaging analysis of system o and use equivariant 
normal form theory in the regime uj = 0(1) to demonstrate that there exists a smooth function a = 
ttTBib, oj, e) at which the system possesses a torus bifurcation in which a stable two-torus is born. First, 
we change variables so that the fold is located at the origin: (x, y) = (1 — x, y -|- |) and a = a — 1. The 
forced vdP equation transforms to 

~/ ~ ~2 ^ 

X = —y + X -X 

" 3 

y' = e{x + a + h cos 9) 

9’ = u. 

Then, to carry out the second-order averaging analysis, it is natural to use the following scaling, which 
comes from the central chart of the desingularization, or blow-up, analysis used in Sectionj^and in |l36l: 

X = \/ex, y = ey, b = \/eb, a = y/ea, 


and to rescale time by t —)■ oot. Hence, after dropping the bars, the system has the following non- 
autonomous form: 


X = - {-y + X ) - —x'^ 
io Su 

y = — {x + a + b cos(to + t)) 

UJ 


(71) 


where S = y/s. The choice of to has no effect on the analysis. 

Next, we apply the near-identity change of variables used in second-order averaging Il4^ . so that 
system ([7T|) transforms into 


6 = — (-6 + ^ 1 ) - b ) 

UJ oUJ 

^2 = — ('Cl + o) + 0{S^), 

UJ 


with G(Ci, C 2 , t, 5) = 0{5^). We label (Tl) as the ‘intermediate’ system; it is smoothly conjugate to the 
original system. 

We are interested in the averaged system. 


X = - i-y + x^) - —X 
(jj 3 lo 

■ ^ / 

y = - {x +a). 

UJ 


;^3 


(73) 


This system has a unique equivariant normal form ll25l due to its symmetry properties. Moreover, 
the time-T map of this normal form must be the normal form of the time-T map, since the two 
operations commute and since equivariant normal forms are unique. Let s denote the Poincare map 
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of this normal form. At a = 0, the eigenvalues of the map s satisfy the non-resonance condition: they are 
not equal to the first four strong resonant eigenvalues. Also, the second Liapunov coefficient is negative; 
in fact, for the averaged system, the second Liapunov coefficient is known to be K6, where iiT < 0 is a 
constant. Hence, at a = 0, the map s satisfies fhe basic hypofheses of fhe Hopf bifurcafion for maps; see 
conditions A) and B), respecfively, and Theorem 2 of |[38l . Therefore, af a = 0, fhe map s undergoes a 
non-degenerafe, supercrifical Hopf bifurcafion in which a normally-hyperbolic invarianf circle is creafed. 
Hence, one also sees direcfly fhaf fhe averaged sysfem ( [7^ undergoes a forus bifurcafion af a = 0 in 
which fhe limif cycle becomes unsfable and a sfable invarianf forus is creafed. 

We now demonsfrafe fhaf fhe full sysfem ( fTT] ) undergoes a forus bifurcafion af some cttb near a = 0, 
in which a sfable invarianf forus is creafed. In parficular, we show fhaf fhe Poincare map s of fhe full 
sysfem ( [TT] ) possesses an invarianf circle 5-close fo fhe one of s. 

Firsf, we observe fhaf fhe same near-idenfify coordinafe change employed above fo puf fhe averaged 
sysfem (731 info ifs equivarianf normal form also pufs fhe infermediafe sysfem (72i info ifs equiv- 
arianf normal form, up fo and including 0(5^). Lef 5 i{z) denofe fhe Poincare map of fhis normal form. 
Again, fhis map musf be fhe fime-T map of fhe normal form, due fo uniqueness in fhe equivarianf 
case. Nexf, we observe fhaf, by standard second-order averaging theory, the Poincare maps are close, 
i.e., on time intervals of length 0 (1/6), 


\5i{z) -s(z)| = 0{5). 

In fact, for each a sufficiently small, the map s* is part of exactly the type of one-parameter family of 
maps studied in ll38]l . with the map of the averaged system being the ‘unperturbed’ map. Hence, for some 
a near a = 0, the map Sj also undergoes a non-degenerate, super-critical Hopf bifurcation in which a 
normally-attracting invariant circle is created. 

Finally, we observe that the time-T maps s and Sj are smoothly conjugate. Hence, it follows that 
the map s of the original system ( f/T] ) also has a non-degenerate Hopf bifurcation, and therefore that the 
original vector field ( fTT] ) has a forus bifurcafion af some otb near zero, in which an affracfing invarianf 
forus is creafed. This concludes fhe demonsfrafion. 

Remark 7. One boundary of fhe paramefer a for fhe exisfence of an invarianf forus for fhe full sysfem 
is given by fhe birfh of fhe canard regime. 

A.2 Asymptotic Expansion of otb ip, u, e) 

In fhis section, we presenf fhe asymptotic expansion of a = aTBp, co, e) for small b. The unforced vdP 
equation has an equilibrium af (x, y) = {a, /(a)), which undergoes a Hopf bifurcafion af o = 1. This 
corresponds fo a periodic orbif of o af 6 = 0 which undergoes a forus bifurcafion af a = 1. We seek 
periodic solutions of ([T]) as an asymptotic series in b, i.e., lef 

OO 

x{t) = ^b’^Xk{t), y{t) 
k=0 

Subsfifufion info Q yields 

xo = yo- fixo), 
yo = e{a - xq), 

af leading order, which is fhe planar van der Pol equation. This has an equilibrium af (x, y) = (a, /(o)). 
The 0{b^) sysfem is 

Xi=yi- /'(xo)xi, 
yi = e(—xi -|- cosut). 

The 0{b^) sysfem is linear wifh solutions fhaf are linear combinations of coswf, sintuf, and exp(Af), 
where 2A = 1 — ± -^(1 — — 4 e. 


k=Q 
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Remark 8. Note that when a = 1 and w = -y/e, there is a resonance. 


When a damped oscillator is driven with a periodic forcing function, the result may be a periodic 
response at the same frequency as the forcing function (see Figure [T3|^a)). Since the unforced oscillation 
is dissipated due to the damping, it is absent from the steady state behaviour. Thus, we seek periodic 
solutions of the 0{b^) system of period T = ^. The solutions are 

— l) euj sm{tuj) + £ (e — uP) cos(ta;) 

{a? - if + {£ - u}‘^f ’ 

e ((a^ “ l) ^ cos(fa;) + uj (a^ — 2a? — e + + 1) sin(fa;)) 

~ 2 I t 2^(2 ■ 

(a^ — 1) uj^ + [£ — uj^) 

We now compute the stability of a periodic solution of O- Let IJ'y') denote the periodic solution 
and let (u, n) = (x — Xy, y — y^) be a perturbation of this orbit. Then the perturbations evolve according 
to 

where the Jacobian evaluated along (xy, yf is the T-periodic matrix: 

The Floquet multipliers pi and p 2 satisfy 




xi(t) = 


piP 2 = exp 


tr Df{x^,yf dt 


A Neimark-Sacker bifurcation occurs when the multipliers are p = for some p. That is, a torus 
bifurcation occurs when ff tr U/(xy, yf dt = 0. This gives the following relation between a, b, e, and 
u! for the location of the torus bifurcation: 


1 — a 


2 


1 6 ^ 

2 (a2 - l)2a;2 + (e _^2)2 


(74) 


Comparisons between the theoretical prediction above and the results of numerical continuation simula¬ 
tions show good agreement for the indicated parameter regions. Divergence from the above perturbation 
analysis occurs precisely at the resonances (figures not shown). 


Remark 9. Note that ( [741 ) is accurate up to terms of order 0{b^). That is, there are no order 0{b‘^) 
corrections in from the b‘^X 2 terms in the asymptotic expansion due to symmetry. More precisely, 
symmetry considerations give 


—2ab^ / X 2 dt = 0. 


7o 


References 

[1] S. Baer, and T. Emeux, Singular Hopf bifurcation to relaxation oscillations, SIAM J. Appl. Dyn. Syst., 46 
(1986), 721-739. 

[2] G. N. Benes, A. M. Barry, T. J. Kaper, M. A. Kramer, and J. Burke, An elementary model of torus canards, 
Chaos, 21 (2011)023131. 

[3] E. Benoit, J.-L. Callot, E. Diener, and M. Diener, Chasse au canard. Collectanea Mathematicae, 31-32 
(1981), 37-119. 

[4] E. Benoit, Canards et enlacements, Inst. Haul. Etud. Sci. Publ. Math., 72 (1990), 63-91. 


34 





[5] K. Bold, C. Edwards, J. Guckenheimer, S. Guharay, K. Hoffman, J. Hubbard, R. Oliva and W. Weckesser, 
The forced van der Pol equation II: canards in the reduced system, SIAM J. Appl. Dyn. Syst., 2 (2003), 
570-608. 

[6] B. Braaksma, Critical phenomena in dynamical systems of van der Pol type, Ph.D. thesis. University of 
Utrecht, 1993. 

[7] M. Br0ns, M. Krupa and M. Wechselberger, Mixed mode oscillations due to the generalized canard phe¬ 
nomenon, in ’’Bifurcation Theory and Spatio-Temporal Pattern Formation”, Fields Institute Communications, 
49 Amer. Math. Soc., Providence, RI, (2006), 39-63. 

[8] J. Burke, M. Desroches, A. M. Barry, T. J. Kaper, and M. A. Kramer, A showcase of torus canards in 
neuronal bursters, J. Math. Neuro., (2012) 2:3 

[9] M. L. Cartwright and J. E. Littlewood, On non-linear differential equations of the second order: I. The 

equation y — k{l — + y = bXk cos(Af + a); k large, Journal of the London Mathematical Society, 20 

(1945) pp. 180-189. 

[10] M.L. Cartwright, Forced oscillations in nonlinear systems Contrib. to theory of nonlinear oscillations, 
Princeton University Press (Study 20) 149-241, 1950. 

[11] S.N. Chow and J.K. Hale. Methods for Bifurcation Theory. Springer-Verlag, 1982. 

[12] A. Delshams and T. M. Seara, An asymptotic expression for the splitting of separatrices of the rapidly forced 
pendulum. Comm. Math. Phys., 150:3 (1992), 443^63. 

[13] M. Desroches, B. Krauskopf, and H. M. Osinga, The geometry of slow manifolds near a folded node, SIAM 
Journal of Applied Dynamical Systems, 7 (2008), 1131-1162. 

[14] M. Desroches, B. Krauskopf, and H. M. Osinga, Numerical continuation of canard orbits in slow-fast 
dynamical systems. Nonlinearity, 23 (2010), 739-765. 

[15] M. Desroches, J. Guckenheimer, C. Kuehn, B. Krauskopf, H. M. Osinga and M. Wechselberger, Mixed-mode 
oscillations with multiple time scales, SIAM Rev., 54 (2012), pp. 211-288. 

[16] M. Desroches, M. Krupa and S. Rodrigues, Inflection, canards and excitability threshold in neuronal models, 
J. Math. Bio., 67 (2013), 989-1017. 

[17] M. Diener, The canard unchained or how fast-slow systems bifurcate. Math. Intelligencer, 6 (1984), 38^9. 

[18] E. J. Doedel, A. R. Champneys, T. F. Fairgrieve, Y. A. Kuznetsov, K. E. Oldeman, R. C. Paffenroth, B. Sanst- 
ede, X. J. Wang and C. Zhang, AUTO-07P: Continuation and bifurcation software for ordinary differential 
equations. Available from: http://cmvl.cs.concordia.ca/ 

[19] F. Dumortier and R. Roussarie, Canard cycles and center manifolds, Mem. Amer. Math. Soc., 577 (1996). 

[20] F. Dumortier, and R. Roussarie, Geometric singular perturbation theory beyond normal hyperbolicity, in 
’’Multiple Time Scales Dynamical Systems” (ed. C. K. R. T. Jones and A. I. Khibnik), IMA Volumes in 
Mathematics and its Applications, Vol. 122, (2001), pp. 29-64. 

[21] W. Eckhaus, Relaxation oscillations including a standard chase on French ducks. Lecture Notes in Math., 
985 (1983) pp. 449^94. 

[22] I. Erchova and D. J. McGonigle, Rhythms of the brain: An examination of mixed mode oscillation approaches 
to the analysis of neurophysiological data. Chaos, 18 (2008), 015115, 14 pp. 

[23] N. Fenichel, Geometric singular perturbation theory for ordinary differential equations. Journal of Differ¬ 
ential Equations, 31 (1979), pp. 53-98. 

[24] J.E. Flaherty, F.C. Hoppensteadt, Frequency entrainment of a forced van der Pol oscillator. Studies in Appl. 
Math., 58 (1978), pp.5-15. 

[25] M. Golubitsky, I. Stewart, and D.G. Schaeffer. Singularities and Groups in Bifurcation Theory. Springer- 
Verlag, 1988. 

[26] J. Guckenheimer and R Holmes, “Nonlinear Oscillations, Dynamical Systems, and Bifurcations of Vector 
Fields”, Springer, 1983. 

[27] J. Guckenheimer, K. Hoffman and W. Weckesser, The forced van der Pol equation I: the slow flow and its 
bifurcations, SIAM J. Appl. Dyn. Syst., 2 (2003), pp. 1-35. 


35 


[28] J. Guckenheimer, Singular hopf bifurcation in systems with two slow variables, SIAM J. Appl. Dyn. Syst., 
7(2008), 1355-1377. 

[29] R. Haiduc Horseshoes in the forced van der Pol system. Nonlinearity, 22, pp. 213-237. 

[30] X. Han and Q. Bi, Slow passage through canard explosion and mixed-mode oscillations in the forced Van 
der Pol’s equation. Nonlinear Dyn., 68 (2012), pp. 275-283. 

[31] E. M. Izhikevich Neural excitability, spiking and bursting. International Journal of Bifurcation and Chaos, 
10(2000), 11711266. 

[32] E. Izhikevich Synchronization of Elliptic Bursters, SIAM Review, 43 (2001), pp. 315-344. 

[33] C. K. R. T. Jones, Geometric singular perturbation theory, in ’’Dynamical Systems” (ed. R. Johnson), 
Lecture Notes in Mathematics, Springer, New York (1995), pp. 44-120. 

[34] M. A. Kramer, R. D. Traub and N. J. Kopell, New dynamics in cerebellar Purkinje cells: torus canards, 
Phys. Rev. Lett., 101 (2008), 068103. 

[35] M. Krupa and P. Szmolyan. Relaxation Oscillation and Canard Explosion, J. Diff. Eq., 174 (2001), 312-368. 

[36] M. Krupa and P. Szmolyan, Extending geometric singular perturbation theory to nonhyperbolic points -fold 
and canard points in two dimensions, SIAM J. Math. Anal., 33 (2001), pp. 286-314. 

[37] M. Krupa and M. Wechselberger, Local analysis near a folded saddle-node singularity. Journal of Differen¬ 
tial Equations, 248 (2010), pp. 2841-2888. 

[38] O.E. Lanford III. Bifurcation of periodic solutions into invariant tori; the work of Ruelle and Takens. In 
Nonlinear problems in the physical sciences and biology. Springer Berlin Heidelberg, pages 159-192, 1973. 

[39] M. Levi, Qualitative analysis of the periodically-forced relaxation oscillations. Memoirs of the AMS, 32:244 
(1981). 

[40] N. Levinson, A second-order differential equation with singular solutions, Ann. Math., 50:1 (1949), pp. 
127-153. 

[41] E. E. Mishchenko and N. Kh. Rozov, Differential equations with small parameters and relaxation oscillations 
(translated from Russian), Plenum, (1980). 

[42] J. Mitry, M. McCarthy, N. Kopell and M. Wechselberger, Excitable Neurons, Firing Threshold Manifolds 
and Canards, J. Math. Neuro., (2013), 3:12. 

[43] B. van der Pol, A theory of the amplitude of free and forced triode vibrations. Radio Review 1 (1920), pp. 
701-710, 754-762. 

[44] B. van der Pol, On relaxation oscillations, Lond., Edinb., Dublin Phil. Mag. and Jl. of Science, Series 7, 2 
(1926), pp. 978-992. 

[45] B. van der Pol, Forced oscillations in a circuit with non-linear resistance (reception with reactive triode), 
Lond., Edinb., Dublin Phil. Mag. and Jl. of Science, Series 7, 3 (1927), pp. 65-80. 

[46] K.-L. Roberts, J. Rubin, and M. Wechselberger, Averaging, folded singularities and torus canards: explain¬ 
ing transitions between bursting and spiking in a coupled neuron model. Preprint. 

[47] H. Rotstein, M. Wechselberger and N. Kopell, Canard induced mixed-mode oscillations in a medial entorhi- 
nal cortex layer 11 stellate cell model, SIAM J. Appl. Dyn. Syst., 7 (2008), 1582-1611. 

[48] J. Rubin and M. Wechselberger, The selection of mixed-mode oscillations in a Hodgkin-Huxley model with 
multiple timescales. Chaos, 18 (2008), 015105, 12 pp. 

[49] J.A. Sanders and E. Verhulst. Averaging Methods in Nonlinear Dynamical Systems. Springer-Verlag, 1982. 

[50] M. Sekikawa, N. Inaba, T. Yoshinaga, and H. Kawakami, Collapse of Duck Solution in a Circuit Driven by 
an Extremely Small Periodic Force, Elec, and Comm, in Japan, Part 3, 88:4 (2005), 199-207. 

[51] P. Szmolyan and M. Wechselberger, Canards in Journal of Differential Equations, 177 (2001), pp. 
419^53. 

[52] P. Szmolyan and M. Wechselberger, Relaxation oscillations in Journal of Differential Equations, 200 
(2004), pp. 69-104. 


36 



[53] W. Teka, J. Tabak, T. Vo, M. Wechselberger and R. Bertram, The dynamics underlying pseudo-plateau 
bursting in a pituitary cell model. Journal of Mathematical Neuroscience, 1 (2011), 12. 

[54] T. Vo and M. Wechselberger, Canards of folded saddle-node type I, In review. 

[55] M. Wechselberger, Existence and bifurcation of canards in in the case of a folded node, SIAM Journal 
of Applied Dynamical Systems, 4 (2005), pp. 101-139. 

[56] M. Wechselberger, A propos de canards (apropos canards). Transactions of the American Mathematical 
Society, 364 (2012), pp. 3289-3309. 

[57] M. Wechselberger, J. Mitry and J. Rinzel, Canard theory and excitability, in ’’Nonautonomous Dynam¬ 
ical Systems in the Life Sciences”, Lecture Notes in Mathematics, Vol. 2102 (Mathematical Biosciences 
Subseries), 2014. 


37 



